Method for fault detection and diagnosis of a rotary machine

ABSTRACT

Modem rotary machine production requires built-in fault detection and diagnoses. The occurrence of faults, e.g. increased friction or loose bonds has to be detected as early as possible. Theses faults generate a nonlinear behavior. Therefore, a method for fault detection and diagnosis of a rotary machine is presented. Based on a rotor system model for the faulty and un-faulty case, subspace-based identification methods are used to compute singular values that are used as features for fault detection. The method is tested on an industrial rotor balancing machine.

CROSS REFERENCE TO RELATED APPLICATIONS

Applicants claim priority under 35 U.S.C. §119 of European ApplicationNo. 07017478.4 filed on Sep. 6, 2007.

BACKGROUND OF THE INVENTION

The invention relates to a method for fault detection and diagnosis of arotary machine, in particular a balancing machine, wherein a rotorhaving an imbalance is rotated and excites a vibration in the rotarymachine due to the imbalance-caused force, and wherein the rotationalspeed of the rotor and the vibrations are measured in order to obtaininput data quantitative for the rotational speed and the vibrations.

SUMMARY OF THE INVENTION

A method for fault detection and diagnosis of a rotary machine, inparticular a balancing machine is provided. A rotor having an imbalanceis rotated and excites a vibration in the rotary machine due to theimbalance-caused force. The rotational speed of the rotor and thevibrations are measured in order to obtain input data quantitative forthe rotational speed and the vibrations.

In accordance with the method, it is assumed that a process of dynamicbehavior of the rotary machine can be modeled by a linear system in theun-faulty state. An over-determined set of linear equations is formed,which contains input and output data of the process and unknown statesof the assumed linear system. The number of states needed to accuratelymodel the process is extracted by using mathematical operations such asorthogonal or oblique projections to form a matrix of which the rankequals the assumed linear system. Singular values are computed by usingSingular Value Decomposition for obtaining an approximate indication forthe order of the assumed linear state.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a rotary system model for normal condition;

FIG. 2 shows the fault state of sliding friction in sensor connection;

FIG. 3 shows the fault state of sliding friction between bearing supportsocket and ground; and

FIG. 4 shows singular values for normal and faulty states.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 1 Introduction

Fault detection and diagnosis is increasingly important for modernrotary machines. Currently, mostly limit checking and periodicmaintenance cycles are used to detect faults. Sometimes signal-basedfault detection is applied. However, these methods mainly rely on theknowledge of experts. This situation can be improved by model-basedfault detection. As more information (proper excitation, process model,several measurements) is used, more accurate fault detection can beperformed. Standard rotary systems behave basically linear. In case ofspecific faults, e.g. sliding friction, loose bonds and motion blocks,linear relations no longer hold. An indication for the nonlinearity canbe used to detect these faults. Taking into account the noisyenvironment the method presented will use subspace approaches toestimate singular values. These singular values can be used to detectthese faults.

2 Modeling

In order to design a model-based fault detection and diagnosis system,the dynamic behavior of the rotor system needs to be modeled. In a firststep, a general model with two degrees of freedom for stiff rotors isgiven. For lower rotary speeds a simplified model can be applied.

It is assumed that the rotor is not fully balanced, so that an imbalanceforce F_(u) and torque M_(u) exist. The rotor is situated on twoindependent bearing supports. Their movement speeds in horizontal planeare denoted by {dot over (x)}₁,{dot over (x)}₂. Plunger coil sensors areused to measure these speeds, resulting in measurement values {dot over(s)}₁,{dot over (s)}₂ (see also FIG. 1).

2.1 Model for Stiff Rotors

Rotor and bearings are assumed to be stiff, the ground connection of thetwo bearing supports is modeled by two spring-damper systems. It isassumed furthermore that the sensor force feedback on the rotor movementcan be neglected.

m_(r)>>m_(s1) m_(r)>>m_(s2)

c₁>>c_(s1) c₂>>c_(s2)

As long as the system stays in normal condition, it can be described bya linear state space system.

{dot over (x)} _(m)(t)=A _(m) x _(m)(t)+B _(m) u _(m)(t)  (1)

y _(m)(t)=C _(m) x _(m)(t)

Applying Newton's law of motion for a rotary mass it follows

$\begin{matrix}{{u_{m}(t)} = \left( {{F_{u}(t)}\mspace{14mu} {M_{u}(t)}} \right)^{T}} & (2) \\{{x_{m}(t)} = \left( {{{\overset{.}{x}}_{1}(t)}\mspace{14mu} {x_{1}(t)}\mspace{11mu} {\overset{.}{x}\;}_{2}(t)\mspace{14mu} {x_{2}(t)}} \right)^{T}} & (3) \\{{y_{m}(t)} = \left( {{{\overset{.}{s}}_{1}(t)}\mspace{14mu} {{\overset{.}{s}}_{2}(t)}} \right)^{T}} & (4) \\{A_{m} = {\frac{1}{m_{r}\theta_{r}}\begin{pmatrix}\left( {{{- \theta_{r}}d_{1}} - {m_{r}d_{1}l_{1}^{2}}} \right) & \left( {{{- m_{r}}c_{1}l_{1}^{2}} - {\theta_{r}d_{1}}} \right) \\0 & {m_{r}\theta_{r}} \\\left( {{l_{2}m_{r}l_{1}d_{1}} - {\theta_{r}d_{1}}} \right) & \left( {{l_{2}m_{r}l_{1}c_{1}} - {\theta_{r}c_{1}}} \right) \\0 & 0 \\\left( {{{- \theta_{r}}d_{2}} + {m_{r}l_{2}d_{2}l_{1}}} \right) & \left( {{{- \theta_{r}}c_{2}} + {m_{r}l_{2}c_{2}l_{1}}} \right) \\0 & 0 \\\left( {{{- l_{2}^{2}}m_{r}d_{2}} - {\theta_{r}d_{2}}} \right) & \left( {{{- \theta_{r}}c_{2}} - {l_{2}^{2}m_{r}c_{2}}} \right) \\0 & {m_{r}\theta_{r}}\end{pmatrix}}} & (5) \\{B_{m} = \begin{pmatrix}{1/m_{r}} & {{- l_{1}}/\theta_{r}} \\0 & 0 \\{1/m_{r}} & {l_{2}/\theta_{r}} \\0 & 0\end{pmatrix}} & (6) \\{C_{m} = \begin{pmatrix}1 & 0 & 0 & 0 \\0 & 0 & 1 & 0\end{pmatrix}} & (7)\end{matrix}$

However, this detailed model is not needed for the regarded rotarysystem. As it is operated in sub-critical region, certainsimplifications can be made.

2.2 Reduced Model

Under the assumption that the machine is driven with sub-critical rotaryspeed, i.e. (ω_(r)<<ω_(crit)) where ω_(r) is the actual rotary speed andω_(crit) is the lower of the two critical speeds according to system (1)it can be assumed that

m₁{umlaut over (x)}₁<<d₁{dot over (x)}₁<<c₁x₁

m₂{umlaut over (x)}₂<<d₂{dot over (x)}₂<<c₂x₂

With these simplifications the model reduces to a model of order two:

$\begin{matrix}{{{\overset{.}{x}}_{r}(t)} = {{A_{r}{x_{r}(t)}} + {B_{r}{u_{r}(t)}}}} & (8) \\{{y_{r}(t)} = {C_{r}{x_{r}(t)}}} & (9) \\{{u_{r}(t)} = \left( {{{\overset{.}{s}}_{1}(t)}\mspace{14mu} {{\overset{.}{s}}_{2}(t)}} \right)^{T}} & (10) \\{{x_{r}(t)} = \left( {{x_{1}(t)}\mspace{14mu} {x_{2}(t)}} \right)^{T}} & (11) \\{{y_{r}(t)} = \left( {{F_{u}(t)}\mspace{14mu} {M_{u}(t)}} \right)^{T}} & (12) \\{A_{r} = \begin{pmatrix}0 & 0 \\0 & 0\end{pmatrix}} & (13) \\{B_{r} = \begin{pmatrix}1 & 0 \\0 & 1\end{pmatrix}} & (14) \\{C_{r} = \begin{pmatrix}c_{1} & c_{2} \\{{- c_{1}}l_{1}} & {c_{2}l_{2}}\end{pmatrix}} & (15)\end{matrix}$

The system is observable and controllable, the poles are on thestability limit. Inputs and outputs are exchanged in order to match thestate space structure. Discretization with small sampling time T₀ leadsto

{dot over (x)} _(d)(k+1)=A _(d) x _(d)(k)+B _(d) u _(d)(k)  (16)

y _(d)(k)=C _(d) x _(d)(k)  (17)

$\begin{matrix}{{u_{d}(k)} = \left( {{{\overset{.}{s}}_{1}(k)}\mspace{14mu} {{\overset{.}{s}}_{2}(k)}} \right)^{T}} & (18) \\{{x_{d}(k)} = \left( {{x_{1}(k)}\mspace{14mu} {x_{2}(k)}} \right)^{T}} & (19) \\{{y_{d}(k)} = \left( {{F_{u}(k)}\mspace{14mu} {M_{u}(k)}} \right)^{T}} & (20) \\{A_{d} = \begin{pmatrix}1 & 0 \\0 & 1\end{pmatrix}} & (21) \\{B_{d} = \begin{pmatrix}T_{0} & 0 \\0 & T_{0}\end{pmatrix}} & (22) \\{C_{d} = \begin{pmatrix}c_{1} & c_{2} \\{{- c_{1}}l_{1}} & {c_{2}l_{2}}\end{pmatrix}} & (23)\end{matrix}$

This model will be used as basis for fault detection.

2.3 Nonlinear Fault States

The presented method is used to detect two fault states where slidingfriction occurs.

2.3.1 Sliding Friction in Sensor Connection

If the sensor is not connected properly to the left moving rotorsupport, the force is propagated via sliding friction. The propagatedforce F_(Rs)=f({dot over (x)}₁−{dot over (s)}₁) is modeled as Coulombdry friction. The sensor dynamics are described by the frequencyresponse G_(s)(w).

$\begin{matrix}{F_{Rs} = {F_{cs}{{sign}\left( {{\overset{.}{x}}_{1} - {\overset{.}{s}}_{1}} \right)}}} & (24) \\{{G_{s}(w)} = \frac{s}{{m_{s}s^{2}} + c_{s}}} & (25)\end{matrix}$

The sensor force feedback on the rotor movement is neglected. Thus, therelation is nonlinear in the outputs equation only.

$\begin{matrix}{{{{\overset{.}{x}}_{p}(t)} = {{A_{p}{x_{p}(t)}} + {B_{p}{u_{p}(t)}}}}{{y_{p}(t)} = \begin{pmatrix}{f\left( {{{\overset{.}{x}}_{1}(t)},{G_{s}(w)},F_{cs}} \right)} \\{{\overset{.}{x}}_{2}(t)}\end{pmatrix}}} & (26) \\{{u_{p}(t)} = \left( {{F_{u}(t)}\mspace{14mu} {M_{u}(t)}} \right)^{T}} & (27) \\{{x_{p}(t)} = \left( {{{\overset{.}{x}}_{1}(t)}\mspace{14mu} {x_{1}(t)}\mspace{14mu} {{\overset{.}{x}}_{2}(t)}\mspace{14mu} {x_{2}(t)}} \right)^{T}} & (28) \\{{y_{p}(t)} = \left( {{{\overset{.}{s}}_{1}(t)}\mspace{14mu} {{\overset{.}{s}}_{2}(t)}} \right)^{T}} & (29) \\{A_{p} = A_{m}} & (30) \\{B_{p} = B_{m}} & (31)\end{matrix}$

2.3.2 Sliding Friction in Ground Connection

If the left rotor support is not properly connected to the ground, thebearing support socket (mass m_(g)) may move on the ground. FIG. (3)shows the dynamic behavior of this fault state. It is assumed that dryCoulomb friction persists between rotor support socket and the ground.

The dynamics may be described by a linear system with nonlinear feedbackaccording to FIG. (3).

{dot over (x)} ₁(t)=A ₁ x ₁(t)+B ₁ u ₁(t)

y ₁(t)=C ₁ x ₁(t)  (32)

with

$\begin{matrix}{\mspace{79mu} {{u_{1}(t)} = \left( {{F_{u}(t)}\mspace{14mu} {M_{u}(t)}\mspace{14mu} {F_{R\; 1}(t)}} \right)^{T}}} & (33) \\{\mspace{79mu} {{x_{1}(t)} = \left( {{{\overset{.}{x}}_{1}(t)}\mspace{14mu} {x_{1}(t)}\mspace{14mu} {{\overset{.}{x}}_{2}(t)}\mspace{14mu} {x_{2}(t)}\mspace{14mu} {{\overset{.}{x}}_{g}(t)}\mspace{14mu} {x_{g}(t)}} \right)^{T}}} & (34) \\{\mspace{79mu} {{y_{1}(t)} = \left( {{{\overset{.}{s}}_{1}(t)}\mspace{14mu} {{\overset{.}{s}}_{2}(t)}} \right)^{T}}} & (35) \\{A_{1} = \left( \begin{matrix}\; & \; & \; & \; & \frac{\left( {{m_{r}d_{1}l_{1}^{2}} + {\ominus_{r}d_{1}}} \right)}{m_{r} \ominus_{r}} & \frac{\left( {{m_{r}c_{1}l_{1}^{2}} + {\ominus_{r}c_{1}}} \right)}{m_{r} \ominus_{r}} \\\; & A_{m} & \; & \; & 0 & 0 \\\; & \; & \; & \; & \frac{\left( {{{- l_{2}}m_{r}l_{1}d_{1}} + {\ominus_{r}d_{1}}} \right)}{m_{r} \ominus_{r}} & \frac{\left( {{{- l_{2}}m_{r}l_{1}c_{1}} + {\ominus_{r}c_{1}}} \right)}{m_{r} \ominus_{r}} \\\; & \; & \; & \; & 0 & 0 \\\frac{c_{1}}{m_{g}} & \frac{d_{1}}{m_{g}} & 0 & 0 & \frac{- c_{1}}{m_{g}} & \frac{- d_{1}}{m_{g}} \\0 & 0 & 0 & 0 & 1 & 0\end{matrix} \right)} & (36) \\{\mspace{79mu} {B_{1} = \begin{pmatrix}\; & \; & 0 \\B_{m} & \; & 0 \\\; & \; & 0 \\\; & \; & 0 \\0 & 0 & {1/m_{g}} \\0 & 0 & 0\end{pmatrix}}} & (37) \\{\mspace{79mu} {C_{1} = \left( {C_{m}\begin{matrix}0 & 0 \\0 & 0\end{matrix}} \right)}} & (38) \\{\mspace{79mu} {F_{Ri} = {F_{d}{{sign}\left( {{\overset{.}{x}}_{g}(t)} \right)}}}} & (39)\end{matrix}$

2.3.3 Summary

The two described fault states introduce nonlinear behavior into thestate space relation, either in the output equation or in the systemequation. An approximation of this behavior with the reduced modelaccording to section 2.2 turns out to be inaccurate.

2.4 Imbalance Force and Torque

It is assumed that the rotor is not fully balanced. The remainingimbalances cause an imbalance force F_(u) and torque M_(u). It isassumed that the rotor speed ω_(r) and rotor roll angle φ_(r) are knownand the imbalance amplitudes A_(u1), A_(u2) and angles φ_(u1), φ_(u2)are measured or known. The imbalance force and torque can be modeled as

$\begin{matrix}{{F_{u}(t)} = {A_{u\; 1}{\omega_{r}^{2}(t)}{\cos \left( {{\phi_{r}(t)} + \phi_{u\; 1}} \right)}}} & (40) \\{\mspace{59mu} {{+ A_{u\; 2}}{\omega_{r}^{2}(t)}{\cos \left( {{\phi_{r}(t)} + \phi_{u\; 2}} \right)}}} & (41) \\{{M_{u}(t)} = {l_{1}A_{u\; 1}{\omega_{r}^{2}(t)}{\cos \left( {{\phi_{r}(t)} + \phi_{u\; 1}} \right)}}} & (42) \\{\mspace{65mu} {{- l_{2}}A_{u\; 2}{\omega_{r}^{2}(t)}{\cos \left( {{\phi_{r}(t)} + \phi_{u\; 2}} \right)}}} & (43)\end{matrix}$

3 Estimation of the Degree of Linearity

As an approximate indication for the degree of linearity, commonsubspace-based methods are well-suited. The mathematical approach thatis used can be described as follows:

-   -   It is assumed that the process can be modeled by a linear state        space system in the un-faulty state.    -   Given the input and output data of the process, an        over-determined set of linear equations is constructed. The set        of equations contains input data, output data and unknown states        of the state space system.    -   The number of states needed to accurately model the process is        extracted by mathematic operations such as orthogonal or oblique        projections. The number of observable and controllable states        equals the system order.    -   The problem of determining what number of states is needed to        model the system is transformed to a matrix rank determination.        This determination is performed approximately by computing the        singular values via Singular Value Decomposition (SVD).

The computed singular values give an approximate indication for theorder of the presumed state space system. If the given process stronglyobeys the reduced model equations according to section (2.2), the methodindicates a process of order two. If nonlinear behavior resides and thelinear model does not fit, the indication becomes indistinct.

3.1 Linearity Indication

This subsection briefly describes the computation of the features forlinearity indication. The used algorithm is partially known as MOESP(Multivariable Output-Error State Space).

3.1.1 State Space System

The input/output relation is assumed to match a linear state spacerelation according to equation (44). N samples of inputs and outputs areavailable.

x(k+1)=Ax(k)+Bu(k)+Bn(k)  (44)

y(k)=Cx(k)+Cm(k)  (45)

The system is observable and controllable of order n. m(k) and n(k)represent white noise sequences. Matrices A, B, C, D and the states x(k)can be transformed by a regular transformation to

Ā=T⁻¹AT

B=T⁻¹B

C=CT

x (k)=T ⁻¹ x(k)  (46)

3.1.2 Alignment in Block Hankel Matrices

The measured data is aligned in block Hankel Matrices

$\begin{matrix}{U_{p} = \begin{pmatrix}{u(0)} & {u(1)} & \ldots & {u\left( {j - } \right)} \\{u(1)} & {u(2)} & \ldots & {u(j)} \\\vdots & \vdots & \ddots & \vdots \\{u\left( { - 1} \right)} & {u()} & \ldots & {u\left( { + j - 2} \right)}\end{pmatrix}} & (47) \\{U_{f\;} = \begin{pmatrix}{u()} & {u\left( { + 1} \right)} & \ldots & {u\left( { + j - 1} \right)} \\{u\left( { + 1} \right)} & {u\left( { + 2} \right)} & \ldots & {u\left( { + j} \right)} \\\vdots & \vdots & \ddots & \vdots \\{u\left( {{2\; } - 1} \right)} & {u\left( {2\; } \right)} & \ldots & {u\left( {{2\; } + j - 2} \right)}\end{pmatrix}} & (48) \\{Y_{p} = \begin{pmatrix}{y(0)} & {y(1)} & \ldots & {y\left( {j - } \right)} \\{y(1)} & {y(2)} & \ldots & {y(j)} \\\vdots & \vdots & \ddots & \vdots \\{y\left( { - 1} \right)} & {y()} & \ldots & {y\left( { + j - 2} \right)}\end{pmatrix}} & (49) \\{Y_{f} = \begin{pmatrix}{y()} & {y\left( { + 1} \right)} & \ldots & {y\left( { + j - 1} \right)} \\{y\left( { + 1} \right)} & {y\left( { + 2} \right)} & \ldots & {y\left( { + j} \right)} \\\vdots & \vdots & \ddots & \vdots \\{y\left( {{2\; } - 1} \right)} & {y\left( {2\; } \right)} & \ldots & {y\left( {{2\; } + j - 2} \right)}\end{pmatrix}} & (50)\end{matrix}$

with

u(k)=({dot over (s)} _(i)(k){dot over (s)} ₂(k))^(T)

y(k)=(F _(u)(k)M _(u)(k))^(T)

N number of measurements 2i maximum order that can be indicated.User-chosen.

j=N−2i+1 if all measurements are used.

The matrices contain all available data and therefore all availableinformation. A set of linear equations is formed which contains theseHankel Matrices and the state Vectors x(k). To explain the procedure,the noise influence is set to zero at this stage.

Y _(p)=Γ_(i) X _(p) +H _(i) U _(p)

Y _(f)=Γ_(i) X _(f) +H _(i) U _(f)  (51)

X _(f) =A ^(i) X _(p)+Δ_(i) U _(p)

To develop this set of equations, following matrices are used:

$\begin{matrix}{\Gamma_{i} = \left( {C^{T}\mspace{14mu} ({CA})^{T}\mspace{14mu} \left( {CA}^{2} \right)^{T}\mspace{14mu} \ldots \mspace{14mu} \left( {CA}^{i - 1} \right)^{T}} \right)^{T}} & (52) \\{\Delta_{i}\left( {A^{i - 1}B\mspace{14mu} A^{i - 2}B\mspace{14mu} \ldots \mspace{14mu} {AB}\mspace{14mu} B} \right)} & (53) \\{H_{i} = \begin{pmatrix}D & 0 & 0 & \ldots & 0 \\{CB} & D & 0 & \ldots & 0 \\{CAB} & {CB} & D & \ldots & 0 \\\vdots & \vdots & \vdots & \ldots & \ldots \\{{CA}^{i - 2}B} & {{CA}^{i - 3}B} & {{CA}^{i - 4}B} & \ldots & D\end{pmatrix}} & (54)\end{matrix}$

The state matrices X_(p) and X_(f) are defined analogously to theinput/output Hankel Matrices:

X _(p)=(x(0)x(1) . . . x(j−1))  (55)

X _(f)=(x(i)x(i+1) . . . x(i+j−1))  (56)

The set of equations (51) can easily be verified by direct insertion. Byremoving the unknown states from this set of equations the solution forY_(f) yields:

$\begin{matrix}{Y_{f} = {\left( {H_{f}\; {\Gamma_{i}\left( {\Delta_{i} - {A^{i}\Gamma_{i}^{\dagger}H_{i}}} \right)}\Gamma_{i}A^{i}\Gamma_{i}^{\dagger}} \right)\begin{pmatrix}U_{f} \\U_{p} \\Y_{p}\end{pmatrix}}} & (57)\end{matrix}$

The notation † stands for the Moore-Penrose-Pseudoinverse. Equation (57)yields direct information on the linearity indication features. If themodel is purely linear, the row space of Y_(f) can be fully described bythe row spaces

$\begin{pmatrix}U_{f} \\U_{p} \\Y_{p}\end{pmatrix}.$

The row space of a matrix is the space spanned by its row vectors. If amatrix J is of full rank, its row space equals the row space of J_(z) ifJ=J_(s)J_(z).

3.1.3 Feature Extraction

For order extraction many different methods are known. The most commonare N4SID, MOESP and CVA. As the underlying system is on the stabilitylimit, the algorithm with the most direct order computation, MOESP, isused. Tests with real data as described in the following have approvedthis choice. MOESP uses a direct RQ decomposition of aligned BlockHankel Matrices:

$\begin{matrix}{{\begin{pmatrix}U_{f} \\U_{p} \\Y_{p}\end{pmatrix} = {\begin{pmatrix}R_{U_{f}} & 0 & 0 & 0 \\0 & R_{U_{p}} & 0 & 0 \\0 & 0 & R_{Y_{p}} & 0 \\{H_{f}R_{U_{f}}} & {{\Gamma_{i}\left( {\Delta_{i} - {A^{i}\Gamma_{i}^{\dagger}H_{i}}} \right)}R_{U_{p}}} & {\Gamma_{i}A^{i}\Gamma_{i}^{\dagger}R_{Y_{p}}} & 0\end{pmatrix}\begin{pmatrix}Q_{U_{r}} \\Q_{U_{p}} \\Q_{Y_{p}} \\Q_{Y_{f}}\end{pmatrix}}}\mspace{79mu} {with}\mspace{14mu} \mspace{79mu} {{U_{f} = {R_{U_{f}}Q_{U_{f}}}},{U_{p} = {R_{U_{p}}Q_{U_{p}}}},{Y_{p} = {R_{Y_{p}}Q_{Y_{p}}}}}} & (58)\end{matrix}$

From the first part, a matrix β_(r), with rank=system order isextracted.

$\begin{matrix}\begin{matrix}{\beta_{r} = \left( {{\Gamma_{i}\left( {\Delta_{i} - {A^{i}\Gamma_{i}^{t}H_{i}}} \right)}R_{U_{p}}\mspace{14mu} \Gamma_{i}A^{i}\Gamma_{i}^{\dagger}R_{Y_{p}}} \right)} \\{= {\left( {\Gamma_{i}\mspace{14mu} \Gamma_{i}} \right)\begin{pmatrix}\left( {\Delta_{i} - {A^{i}\Gamma_{i}^{t}H_{i}}} \right) & 0 \\0 & {A^{i}\Gamma_{i}^{t}R_{Y_{p}}}\end{pmatrix}}}\end{matrix} & \begin{matrix}(59) \\(60)\end{matrix}\end{matrix}$

The rank of β_(r) equals the number of linear independent vectors in itscolumn space, which means that

rank(β_(r))=rank(Γ_(i))=n  (61)

The column space of a matrix is the space spanned by its column vectors.If a matrix J is of full rank, its column space equals the column spaceof J_(z) if J=J_(s)J_(z).

A proper method can extract the ‘true’ order of the underlying system.In the case of a faulty, nonlinear behavior, the linearity indicationdiffers from the described model of order two. To extract the matrixrank in an approximate way, Singular Value Decomposition (SVD) is used.The SVD of β_(r) yields 3 matrices U₁, S₁, V₁

$\begin{matrix}\begin{matrix}{\beta_{r} = {U_{1}S_{1}V_{1}^{T}}} \\{= {{U_{1}\begin{pmatrix}\sigma_{1} & 0 & \ldots & 0 \\0 & \sigma_{2} & \ldots & 0 \\\vdots & \vdots & \ddots & \vdots \\0 & 0 & \ldots & \sigma_{i}\end{pmatrix}}V_{1}^{T}}}\end{matrix} & \begin{matrix}(62) \\(63)\end{matrix}\end{matrix}$

U, and V₁ are orthogonal matrices. S₁ is a diagonal matrix whichcontains the singular values σ_(i). In case of a fault-free, not noisysystem, n singular values are nonzero while all other singular valuesare zero.

σ₁> . . . >σ_(n)>0  (64)

σ_(n+1)= . . . =σ_(2i)=0  (65)

3.2 Noise Influence on Feature Extraction

Under the influence of noise (process noise as well as measurementnoise), the SVD no longer yields clear order decisions. To representnoise influences, it will be assumed that the measurements y(k) arecontaminated by white noise n(k)=(n₁(k)n₂(k))^(T). The excitations u(k)contain noise m(k)=(m₁(k)m₂(k))^(T). With

$N_{p} = \begin{pmatrix}{n(0)} & {n(1)} & \ldots & {n\left( {j - } \right)} \\{n(1)} & {n(2)} & \ldots & {n(j)} \\\vdots & \vdots & \ddots & \vdots \\{n\left( { - 1} \right)} & {n()} & \ldots & {n\left( { + j - 2} \right)}\end{pmatrix}$ $N_{f} = \begin{pmatrix}{n()} & {n\left( { + 1} \right)} & \ldots & {n\left( { + j - 1} \right)} \\{n\left( { + 1} \right)} & {n\left( { + 2} \right)} & \ldots & {n\left( { + j} \right)} \\\vdots & \vdots & \ddots & \vdots \\{n\left( {{2\; } - 1} \right)} & {n\left( {2\; } \right)} & \ldots & {n\left( {{2\; } + j - 2} \right)}\end{pmatrix}$

and subsequent formation of M_(p) and M_(f), equation (51) is enhancedto

Y _(p)=Γ_(i) X _(p) +H _(i) U _(p) +H _(i) N _(p) +M _(p)

Y _(f)=Γ_(i) X _(f) +H _(i) U _(f) +H _(i) N _(f) +M _(f)  (66)

X _(f) =A ^(i) X _(p)+Δ_(i) U _(p)+Δ_(i) N _(p)

The solution for Y_(f) then yields

$\begin{matrix}{Y_{f} = {\left( {H_{f}\mspace{14mu} {\Gamma_{i}\left( {\Delta_{i} - {A_{\mspace{14mu}}^{i}\Gamma_{i}^{\dagger}H_{i}}} \right)}\; \Gamma_{i}A^{i}\Gamma_{i}^{\dagger}} \right)\begin{pmatrix}U_{f} \\U_{p} \\Y_{p}\end{pmatrix}}} & (67) \\{{+ \left( {{H_{i}\mspace{14mu} I\mspace{14mu} {\Gamma_{i}\left( {\Delta_{i} - {A_{\mspace{14mu}}^{i}\Gamma_{i}^{\dagger}H_{i}}} \right)}} - {\Gamma_{i}A^{i}\Gamma_{i}^{\dagger}}} \right)}\begin{pmatrix}N_{f} \\M_{f} \\N_{p} \\M_{p}\end{pmatrix}} & (68)\end{matrix}$

and the SVD is performed on

$\begin{matrix}{\beta_{i} = {\beta_{r} + {\left( {{H_{i}\mspace{14mu} I\mspace{20mu} {\Gamma_{i}\left( {\Delta_{i} - {A^{i}\Gamma_{i}^{t}H_{i}}} \right)}} - {\Gamma_{i}A^{i}\Gamma_{i}^{\dagger}}} \right){\begin{pmatrix}N_{f} \\M_{f} \\N_{p} \\M_{p}\end{pmatrix}/U_{f}}}}} \\{= {\beta_{r} + E_{n}}}\end{matrix}\begin{pmatrix}U_{p} \\Y_{p}\end{pmatrix}$

where

β_(i)=2i×2i matrix with rank that has to be estimated

β_(r)=2i×2i matrix with rank n (n=system order)

E _(N)=2i×2i Matrix with full rank (noise representation)

The notation

$/{U_{f}\begin{pmatrix}U_{p} \\Y_{p}\end{pmatrix}}$

stands for the orthogonal projection onto the row space of

$\begin{pmatrix}U_{f} \\U_{p} \\Y_{p}\end{pmatrix}.$

where only the part lying in the row space of

$\quad\begin{pmatrix}U_{p} \\Y_{p}\end{pmatrix}$

is considered. In literature, this projection is referred to as ‘ObliqueProjection’.

It is shown in [VOdM96] that E_(n)→0 for N→∞.

An infinite number of measurements is not achievable. However, if theSignal-to-Noise-Ratio (SNR) and the measurement number N is sufficientlyhigh, the system impact on the Singular Values is larger than the noiseinfluence (see [Lju99]) and E_(n)<<β_(r) holds and therefore rank(β_(i))=rank (β_(r)). Normally, the influence of E_(n) is not fullynegligible and the Singular Values computation yields

σ₁> . . . >σ_(n)>>σ_(n+1), . . . , σ_(2i)  (69)

The singular values representing the system structure dominate, allsuccessive singular values represent the noise influence and areconsiderably lower.

4 Results

As test rig, a industrial production machine for rotors with massm_(r)≈25 kg is used. It is equipped with 2 standard plunger coilsensors. The fundamentation can be adjusted by common clamps. Themachine is driven in sub-critical rotary speed. The environment consistsof normal industrial surrounding, e.g. other machines, noise etc.

Three different runs are analyzed:

-   -   Un-faulty standard run    -   Run with sliding friction in sensor connection    -   Run with sliding friction in ground connection

The features described in the preceding section are computed for the 3described states. All tests are done on the same machine with equalbearings and rotor. For each run, a timespan of 5 seconds is examined.The singular values are computed according to the preceding chapter.FIG. (4) gives an overview over the different features. The values aregiven in logarithmic scale, the standard deviation over all regardedruns is indicated.

The first two Singular Values, which represent the linear behavior,remain nearly equal. The third and fourth singular value refer to modelinaccuracies and noise. In case of nonlinear behavior, their values areconsiderably higher than in the purely linear, un-faulty case.Principally, these values can be used as a feature for the appearance ofnonlinear faults in linear systems.

The invention described above presents a subspace-based method to detectthe occurrence of nonlinear fault states in linear systems. A rotorsystem was used as example. The fault-free state as well as twodifferent fault cases have been modeled and tested on an industrialrotor balancing machine. It has been shown that the computed singularvalues are highly sensible to nonlinear fault states and are well-suitedas features for the occurrence of the considered friction faults.

Symbols and Abbreviations Symbol Unit Description F_(u) [N]imbalance-caused force M_(u) [Nm] imbalance-caused torque {dot over(x)}₁ [m/sec] movement speed of left rotor support {dot over (x)}₂[m/sec] movement speed of right rotor support {dot over (s)}₁, {dot over(s)}₂ [m/sec] measurements by speed sensors c₁, c₂ [N/m] springconstants d₁, d₂ [N/m] damper constants m_(r) [kg] rotor mass Θ_(r)[kgm²] rotor moment of inertia l₁, l₂ [m] distance to the center ofgravity A_(u1), A_(u2) [l] imbalance amplitude φ_(u1), φ_(u2) [l]imbalance amplitude ω_(r) [rad/sec] rotation speed φ_(r) [rad] rotorroll angle n(k), m(k) white noise sequences

FIG. 1: Rotary system model for normal condition. The rotor (mass m_(r))movement in horizontal plane can be modeled with linear spring-dampersystems (left drawing, birdview). The sensor plunger coils (massesm_(s1),m_(s2)) are connected directly to the bearings. As example, theconnection of plunger coil 1 is given in the right drawing (side view).

FIG. 2: Fault state of sliding friction in sensor connection

FIG. 3: Fault state of sliding friction between bearing support socketand ground.

FIG. 4: Singular values for normal and faulty states

1. Method for fault detection and diagnosis of a rotary machine, inparticular a balancing machine, wherein a rotor having an imbalance isrotated and excites a vibration in said rotary machine due to theimbalance-caused force, and wherein the rotational speed of said rotorand said vibrations are measured in order to obtain input dataquantitative for said rotational speed and said vibrations, the methodfurther comprising: assuming that a process of dynamic behavior of saidrotary machine can be modeled by a linear system in the un-faulty state;forming an over-determined set of linear equations which contains inputand output data of the process and unknown states of the assumed linearsystem; extracting the number of states needed to accurately model theprocess by using mathematic operations such as orthogonal or obliqueprojections to form a matrix of which the rank equals the assumed linearsystem; computing singular values by using Singular Value Decompositionfor obtaining an approximate indication for the order of the assumedlinear system.
 2. The method of claim 1, wherein the assumed linearsystem is formed as a linear state space system.
 3. The method of claim1, wherein the common MOESP method which uses a direct decomposition ofaligned Block Hankel Matrices is used for order extraction of the numberof states.
 4. The method of claim 1, wherein the rotor is rotated with asub-critical rotational speed.
 5. Computer program using the method ofclaim 1.